fileR <- dir("../script", ".R$")
for (ii in fileR)
source(paste0("../script/", ii))
data <- read.table("../data/data_mut_DD.csv", header = TRUE, as.is = TRUE)
head(data)
t1 <- as.numeric(data$mut_lof)
t1 <- t1[t1 > 0]
data[data$mut_lof == 0, ]$mut_lof <- min(t1)/10
#data <- data[data$mut_missense >0, ]
t2 <- as.numeric(data$mut_missense)
t2 <- t2[t2 > 0]
data[data$mut_missense == 0, ]$mut_missense <- min(t2)/10
dim(data)
#data <- data[data$mut_missense >0, ]
t2 <- as.numeric(data$mut_damaging)
t2 <- t2[t2 > 0]
data[data$mut_damaging == 0, ]$mut_damaging <- min(t2)/10
allDNData <- data[, paste0("dn_", c("damaging", "lof"), "_DD")]
allMutData <- data[,paste0("mut_", c("damaging", "lof"))]
head(data.frame(allMutData, allDNData))
ntrioDD = 4293
N <- list(dn = ntrioDD)
modelDNdata <- list(NN = dim(allDNData)[1], #Gene numbers
K = 2, #Hypothesis numbers: should be default
NCdn = 2, #Number of de novo classes
Ndn = array(rep(N$dn, 2)), # Family numbers
dataDN = data.frame(allDNData), # Denovo data
mutRate = data.frame(allMutData), # Mutation rates
betaPars = c(6.7771073, -1.7950864, -0.2168248), #Adjust beta's values: should be default
adjustHyperBeta = as.integer(1), ##1 if want to adjust beta, 0 if don't want to adjust beta
upperPi0 = 0.5, lowerPi0 = 0, lowerBeta = 0, ##Lower and upper limits of pi: should be default
lowerHyperGamma = 1, lowerGamma = 1, #Should be default
hyperBetaDN0 = array(c(5, 1)),
hyper2GammaMeanDN = array(c(1, 1)), ##Priors for mean RRs: meanRR ~ gamma(hyper2GammaMeanDN, hyper2BetaDN)
hyper2BetaDN = array(c(0.05, 0.05)) ##Priors for mean RRs: meanRR ~ gamma(hyper2GammaMeanDN, hyper2BetaDN)
)
nIteration = 1000 ##T
nChain = 1
nCore = 1
nThin = floor(nIteration/1000)
mcmcDD <- stan(model_code = DNextTADA, ##Use only De novo model inside extTADA
data = modelDNdata, ##Model data as described above
iter = nIteration, ##
chains = nChain,
cores = nCore,
thin = nThin)
stan_trace(mcmcDD)
mcmcDD
pars0 <- estimatePars(pars = c('pi0',
'hyperGammaMeanDN[1]', 'hyperGammaMeanDN[2]',
'hyperBetaDN[1]', 'hyperBetaDN[2]'),
mcmcResult = mcmcDD)
pars0
par(mfrow = c(1, 2))
plotParHeatmap(pars = c("pi0", "hyperGammaMeanDN[1]"), mcmcResult = mcmcDD)
plotParHeatmap(pars = c("pi0", "hyperGammaMeanDN[2]"), mcmcResult = mcmcDD)
##Get gene list
geneName <- data[, 1]
##Set parameters: use pars0 above
parsFDR <- list(gammaMeanDN = pars0[, 1][2:3],
betaDN = pars0[, 1][4:5],
pi0 = pars0[, 1][1],
nfamily = rep(ntrioDD, 2))
dataFDR <- calculateFDR(pars = parsFDR,
dnData = allDNData, mutData = allMutData,
geneName = geneName)
head(dataFDR)
dim(dataFDR[dataFDR$qvalue < 0.1, ])
dim(dataFDR[dataFDR$qvalue < 0.05, ])